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14.  ABSTRACT 
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(FWUAV).  Flapping  wing  flight  is  a  complex  phenomenon  encompassing  unsteady  effects,  controls  using  multiple  degrees  of 
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phenomena  are  not  well  modeled  using  the  traditional  conceptual  aerodynamic  models  originally  developed  for  fixed  wing 
and  rotary  wing  aircraft.  In  this  study.  Blade  Element  Theory  is  combined  with  momentum  theory  (called  Blade  Element 
Momentum  Theory  [BEMT])  to  estimate  aerodynamic  loads  on  a  FWUAV.  The  BEMT  model  is  also  combined  with 
experimental  scans  of  a  FWUAV  wing  in  a  wind  tunnel  to  represent  the  actual  wing  shape  during  flight  (represented  by  three- 
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this  lift  coefficient  estimation  is  used  in  conjunction  with  BEMT  to  create  a  comprehensive  model  for  flapping  wing  flight  and 
model  calculations  are  compared  against  experimental  data. 
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1.  Introduction/Background 


The  phenomenon  of  flapping  wing  flight  in  nature  has  been  studied  for  centuries.  Recently, 
flapping  flight  for  unmanned  aerial  vehicle  (UAV)  applications  has  become  an  area  of  interest 
for  military  and  civilian  applications  alike.  Flapping  wing  flight  offers  many  potential 
advantages  over  traditional  fixed  and  rotary  wing  aircraft.  Fixed  wing  UAVs  are  a  tried  and  true 
platform  in  theater,  with  platforms  such  as  the  MQ-1  Predator,  the  RQ-11  Raven,  and  the  RQ-7 
Shadow — all  of  which  are  currently  in  use  by  the  U.S.  military.  These  fixed  wing  vehicles  have 
the  advantages  of  long  range  and  endurance,  and  high  payload  capabilities;  however,  they 
require  high  forward  flight  speeds  and  most  configurations  cannot  hover,  which  makes  them 
difficult  to  control  in  confined  spaces. 

Conversely,  rotary  wing  UAVs  are  highly  maneuverable,  can  fly  at  lower  forward  speeds,  and 
can  hover,  but  generally  have  lower  endurance  times  and  are  louder  than  fixed  wing  UAVs  due 
to  high  rotor  tip  speeds.  The  goal  of  using  flapping  wing  UAVs  (FWUAVs)  is  to  bridge  the  gap 
between  fixed  and  rotary  wing  UAVs.  FWUAVs  have  the  potential  to  fly  at  lower  airspeeds  than 
fixed  wing  aircraft  and  most  have  the  ability  to  hover,  which  enables  FWUAVs  to  be  flown  in 
confined  spaces.  Compared  to  rotorcraft,  FWUAVs  tend  to  be  quieter  since  the  wing  flapping 
speed  is  generally  much  slower  than  a  rotor’s  rotation.  This  combination  of  maneuverability, 
hover  capability,  and  stealthiness  makes  FWUAVs  a  potential  choice  for  use  in  confined  spaces. 
Figure  1  shows  an  example  of  each  type  of  UAV. 


Figure  1.  An  example  of  each  types  of  UAV:  (left  to  right)  fixed  wing  UAV  (7),  rotary  wing  UAV  (2), 
and  FWUAV  (7). 

Two  main  configurations  of  FWUAVs  are  currently  being  investigated:  bird-inspired  and  insect- 
inspired  flight.  Bird- inspired  flight  is  based  on  the  forward  posture  of  birds,  where  a  forward 
velocity  is  required  to  maintain  lift  and  the  flapping  motion  is  in  a  roughly  vertical  plane  with 
respect  to  the  forward  velocity.  Conversely,  insect-inspired  flight  is  characterized  by  a  hovering 
posture,  where  insects’  wings  oscillate  in  a  plane  horizontal  to  the  ground.  At  the  beginning  and 
end  of  each  stroke  (flap),  insect  wings  will  rotate,  as  shown  in  figure  2. 
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Figure  2.  Insect  flapping  (5). 

This  project  is  based  on  the  bird-inspired  FWUAV  called  the  RoboRaven  (figure  3),  which  was 
design  and  built  by  researchers  at  the  University  of  Maryland  (UMD)  (4).  It  has  a  wingspan  of 
0.9  m  and  a  max  chord  length  of  ~0.3  m.  The  wings  can  move  independently  of  one  another, 
driven  by  two  separate  servos,  which  rotate  the  wings  at  approximately  3-4  Hz.  In  flight,  the 
RoboRaven  cruises  at  approximately  6.7  m/s,  which  is  in  the  flight  regime  of  Reynolds  number 
-120,000. 


Figure  3.  UMD  RoboRaven. 

The  overall  goal  of  this  project  is  to  improve  upon  the  accuracy  of  existing  FWUAV 
aerodynamic  models  to  be  used  in  the  conceptual  design  process.  FWUAV  design  is  currently  a 
sequential  trial  and  error  process,  where  engineers  iterate  through  many  different  designs  until 
reaching  a  desirable  configuration.  This  process  is  very  labor  and  time  intensive,  since  designers 
must  physically  build  and  test  every  component  on  the  FWUAV.  Aerodynamic  models  that  can 
accurately  predict  flight  forces  enable  designers  to  run  through  many  iterations  of  design  prior  to 
building  any  prototypes.  Additionally,  pairing  the  electric  motor  that  will  move  the  FWUAV ’s 
wings  to  the  flapping  cycle  and  flight  forces  is  very  important  to  efficient  operation  in  an 
FWUAV.  Integrating  these  two  fields  earlier  in  the  design  could  potentially  yield  improved 
designs  and  shorter  design  times. 

In  order  to  accomplish  the  objectives  outlined,  two  main  goals  were  identified: 
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1.  Implement  Blade  Element  Momentum  Theory  (BEMT)  in  the  J.D.  DeLaurier  Blade 
Element  Theory  flapping  wing  aerodynamic  code  (5). 

2.  Use  wing  shape  data  from  Digital  Image  Correlation  (DIC)  tunnel  experiments  of  an  actual 
RoboRaven  UAV  to  calculate  lift  coefficient  across  the  UAV  wing  span: 

a.  Use  thin  airfoil  theory  to  calculate  the  lift  coefficient  (Q) 

The  above-mentioned  BEMT  aerodynamic  analysis  is  created  using  a  combination  of  two 
aerodynamic  theories:  Blade  Element  Theory  (BET)  and  momentum  theory.  In  BET,  the  wing  of 
a  UAV  is  discretized  across  the  entire  span  into  chordwise  “slices.”  Each  of  these  slices 
experiences  the  flight  forces  of  lift,  drag,  and  thrust  (in  differential  form).  Because  the  wing  is 
three  dimensional  (3-D),  each  slice  is  really  a  small  section  with  a  width  in  the  spanwise 
direction;  however,  it  is  treated  as  an  airfoil  with  infinite  span.  Once  all  the  differential  forces 
have  been  calculated  at  each  slice,  integration  is  performed  along  the  entire  span  to  calculate  the 
total  forces  on  the  wing,  as  shown  in  figure  4. 


A I  iff 


IL 


Figure  4.  (Left)  Discretize  wing  (5)  and  wing  slide  forces. 

Compared  to  BET,  momentum  theory  (figure  5)  is  a  much  more  simplified  approach  to 
calculating  aerodynamic  forces.  In  momentum  theory,  the  momentum  change  of  moving  air 
deflected  off  a  wing  is  used  to  calculate  lift  and  thrust.  When  a  wing  is  placed  at  an  angle  of 
attack,  a,  with  respect  to  the  forward  velocity,  U,  it  will  deflect  air  downward  at  a  velocity,  w, 
called  induced  velocity.  Induced  velocity  is  ultimately  used  in  momentum  theory  to  calculate  lift 
and  thrust. 


Figure  5.  Momentum  theory  induced  velocity. 
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For  a  flapping  wing,  momentum  theory  analyzes  the  wings  as  a  “partial  actuator  disk”  (figure  6), 
where  the  disk’s  size  is  determined  by  the  swept  area  of  the  wing  flapping  when  viewing  the 
UAV’s  frontal  area.  The  momentum  change  of  the  air  moving  across  this  partial  disk  area  is  used 
to  calculate  thrust.  Thus,  in  momentum  theory,  the  entire  wing  flapping  motion  is  represented  by 
the  partial  disk  area,  as  opposed  to  discretizing  the  wing  along  the  span  as  in  BET. 


Figure  6.  Momentum  theory  partial  disk  area  (adapted  from 
Shkarayev  and  Silin  [6]). 

In  both  BET  and  momentum  theory,  there  is  one  unknown  variable,  induced  velocity,  in  the 
equations  used  for  thrust  calculation  (see  figure  5).  When  BET  and  momentum  theory  are  used 
separately,  induced  velocity  is  either  approximated  using  a  function  or  is  assigned  a  constant 
value  as  defined  by  the  individual  performing  the  aerodynamic  analysis.  BEMT  solves  this 
problem  by  combining  the  two  equations  for  BET  and  momentum  theory  in  order  to  solve  for 
downwash  velocity.  Once  downwash  velocity  is  calculated,  it  is  used  in  the  BET  analysis  at  each 
spanwise  “slice”  and  forces  are  calculated. 

In  the  DeLaurier  flapping  wing  BET  model  (5),  the  lift  coefficient  is  calculated  using  a  function 
approximation;  however,  throughout  the  flapping  cycle,  a  FWUAV’s  wing  shape  will  change 
dramatically.  It  is  necessary  to  account  for  these  shape  changes  in  an  aerodynamic  modeling 
code,  because  wing  shape  is  extremely  important  when  calculating  the  lift  coefficient.  In  order  to 
account  for  these  shape  changes,  DIC  can  be  used  to  measure  wing  deformation  of  a  FWUAV 
wing  in  a  wind  tunnel.  DIC  uses  high  resolution  pictures  of  a  patterned  wing  surface  to  track 
movements  while  the  wing  is  flapping.  DIC  images  of  the  RoboRaven  (taken  by  UMD 
researchers  [4])  throughout  the  flap  cycle  offer  a  snapshot  of  the  wing  shape,  represented  as  a 
scatter  plot  of  the  wing  (figure  7). 
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Figure  7.  DIC  wing  covered  in  markers  (4). 


2.  Experiment/Calculations 


The  experimental  phase  of  this  project  was  completed  entirely  via  computer  modeling  in 
MATLAB,  thus  the  following  section  describes  the  equations  and  background  for  implementing 
the  modeling  process.  It  is  divided  into  two  sections,  following  the  two  main  goals  outlined  in 
the  background  section:  creation  of  the  blade  element  momentum  theory  model  and  use  of  3-D 
wing  scatter  plots  with  thin  airfoil  theory  to  calculate  the  lift  coefficient. 

2.1  Blade  Element  Momentum  Theory  Model  Creation 

As  described  in  section  1,  the  purpose  of  combining  BET  and  momentum  theory  is  to  calculate 
induced  velocity  using  the  two  theories  together,  as  opposed  to  using  a  functional  approximation 
of  induced  velocity.  J.D.  DeLaurier’s  BET  model  was  used  as  a  basis  for  creating  the  BEMT 
model  in  this  project,  in  which  induced  velocity  is  approximated  using 


Wo  _  2Q0  +  0) 
U  ~  2+AR 

(1) 


(1) 


where  U  =  free  stream  velocity,  a0  —  zero  lift  angle,  9  =  flight  path  angle,  and  AR  =  aspect  ratio. 


In  this  project,  BEMT  is  used  to  improve  upon  DeLaurier’s  induced  velocity  approximation  by 
integrating  the  momentum  model  from  Shkarayev  and  Silin’s  2010  paper  (6)  with  DeLaurier’s 
BET  model.  Shkarayev ’s  momentum  model  estimates  the  thrust  of  an  FWUAV  by  examining  the 
air  moving  through  the  plane  of  wing  flapping,  called  a  “partial  disk  area”  (see  figure  6).  It  is 
assumed  that  the  flow  through  the  partial  disk  area  is  uniform  throughout.  The  partial  disk  area  is 
defined  by  the  flap  angle,  <D>,  and  wing  span,  b,  according  to  equation  2.  One  shortcoming  of  this 
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approach  is  that  the  flow  is  not  uniform  throughout  the  disk,  as  flow  velocity  is  highly  dependent 
on  wing  shape  and  spanwise  position  during  flapping.  A  potential  improved  approach  is 
described  in  section  4,  where  future  work  is  presented. 


A  = 


c Pb2 

~T~ 


(2) 


Using  the  inputs  seen  in  equation  3  and  described  below,  Shkarayev’s  momentum  theory 
equation  will  approximate  either  induced  velocity,  if  thrust  is  known,  or  thrust,  if  the  induced 
velocity  is  input.  Since  it  is  desired  to  find  a  value  for  the  induced  velocity  to  be  used  in  the  BET 
analysis,  thrust  must  be  found  through  another  means.  This  is  accomplished  in  the  BEMT  model 
creation  procedure,  which  is  described  below. 

The  iteration  procedure  consists  of  estimating  thrust,  using  an  initial  “guess”  for  the  thrust  of  the 
FWUAV  (figure  8).  This  thrust  estimate  is  used  to  calculate  induced  velocity  in  equation  3  from 
Shkarayev  and  Silin’s  paper  (6): 
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w4  +  2w3Vcos({?)  +  w2V2  =  ,  (3) 

where  w  =  induced  velocity,  V  =  flight  velocity,  P  =  angle  of  attack,  T  =  thrust,  p  =  air  density, 
and  A  =  wing  partial  disk  area. 


Figure  8.  BEMT  iteration  loop. 

Once  an  initial  calculation  for  induced  velocity  from  momentum  theory  is  obtained,  it  is  used  in 
the  original  DeLaurier  BET  code  in  place  of  the  induced  velocity  function  approximation  shown 
in  equation  3.  The  BET  code  is  used  to  calculate  thrust,  and  this  thrust  is  compared  to  the 
original  thrust  estimate.  Depending  on  the  optimization  routine  used  in  this  iteration  loop,  the 
thrust  estimate  is  systematically  modified  until  the  error  between  the  BET  calculated  thrust  and 
the  estimated  thrust  is  minimized.  This  iteration  loop  should  result  in  agreement  between 
momentum  theory  and  BET,  where  the  induced  velocity  estimate  used  in  momentum  theory  and 
BET  will  result  in  the  smallest  error  between  both  methods. 
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2.2  DIC  Wing  Shape  Measurements 


As  described  previously,  a  FWUAV’s  wing  shape  will  vary  greatly  throughout  the  flap  cycle, 
and  the  lift  coefficient  is  highly  dependent  on  the  wing  shape.  Thus,  it  is  important  to  capture 
wing  shape  changes  during  flapping  in  order  to  have  an  accurate  representation  of  the  lift 
coefficient.  This  project  accomplishes  the  goal  of  calculating  the  lift  coefficient  based  on  wing 
shape  through  the  use  of  DIC  wing  deformation  scatter  plots,  an  example  of  which  is  seen  in 
figure  9.  These  scatter  plots  are  simple  3-D  representations  of  the  wing,  taken  directly  from  DIC 
images  of  an  FWUAV  flapping  in  a  wind  tunnel. 


Figure  9.  (Left)  The  3-D  wing  scatter  plots  throughout  flap  cycle  and  (right)  a  3-D  wing  scatter  plot. 

Once  the  3-D  wing  shape  scatter  plots  have  been  obtained  from  DIC  images,  it  is  necessary  to 
calculate  the  lift  coefficient  to  be  used  in  BET  portion  of  the  BEMT  code.  In  BET,  the  wing  is 
represented  by  “slices”  along  the  span,  thus  it  is  necessary  to  calculate  the  lift  coefficient  at  each 
of  these  spanwise  slices  (figure  10). 


Figure  10.  (Left)  Wing  chord  shapes  along  the  span  and  (right)  two-dimensional  (2-D)  airfoil  shapes  from  the  3-D 
wing  scatter  plots. 

One  method  of  estimating  the  lift  coefficient  is  through  the  use  of  thin  airfoil  theory,  which 
calculates  Q  using  a  polynomial  approximation  of  airfoil  shape.  Q  is  calculated  according  to 
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equation  4  from  Anderson  (7),  where  a  =  angle  of  attack,  dz/dx  =  derivative  of  polynomial 
approximation  of  wing  shape,  and  9  =  chordwise  location  of  airfoil  in  polar  coordinates.  The 
integral  shown  is  in  terms  of  polar  coordinates,  where  one  will  integrate  from  the  leading  edge  of 
the  airfoil  to  the  trailing  edge  of  the  airfoil,  or  0  chord  length  to  1  chord  length.  Since  the 
equation  is  converted  to  polar  coordinates,  the  bounds  become  0  to  jt.  Further  explanation  can  be 
found  in  Anderson  (7). 


c  i  =  2  n\a  +  —  f  —(cos  6  -  l)d  0 ]  (4) 

k  J  dx 
0 

In  order  to  use  thin  airfoil  theory,  it  is  necessary  to  convert  from  3-D  wing  scatter  plots  to 
discrete  airfoils  along  the  span,  where  the  airfoils  are  approximated  by  a  polynomial  function, 
z(x)  (figure  11).  This  was  accomplished  by  grouping  wing  section  slices  from  3-D  scatter  plots 
along  the  span  into  2-D  airfoil  shapes,  as  seen  in  the  chart  in  figure  11.  Depending  on  the  number 
of  slices  (e.g.,  10),  the  wing  would  be  split  into  equal  parts  along  the  span.  Each  of  these  groups 
of  points  is  then  viewed  solely  in  two  dimensions  such  that  the  wing  shape  is  seen  as  a  2-D 
airfoil.  Thus,  each  spanwise  wing  section  is  represented  as  a  simple  2-D  airfoil,  which  is  then 
used  in  conjunction  with  thin  airfoil  theory. 
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Figure  11.  Thin  airfoil  theory  Q  calculation  from  a  3-D  wing  scatter  plot. 


3.  Results  and  Discussion 


3.1  BEMT  Code  Creation 

The  first  phase  of  the  project  was  to  implement  BEMT  using  a  combination  of  the  DeLaurier 
BET  code  and  the  Shkarayev  momentum  theory.  In  the  original  DeLaurier  paper,  a  model 
pterosaur  is  used  as  an  example  for  running  the  code.  The  pterosaur’s  geometric  properties  are 
provided  as  inputs,  such  as  wing  planform  shape  and  span,  as  well  as  flapping  cycle  properties 
such  as  amplitude  of  flap  angle  and  flapping  speed.  In  the  code,  the  amount  of  “(3  =  dynamic 
twist”  (wing  twisting  along  the  span)  was  used  for  analysis,  where  a  more  dynamic  twist  equates 
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to  more  twisting  along  the  wing  span.  As  the  dynamic  twist  was  changed,  the  resulting  lift  and 
thrust  were  analyzed.  The  complete  set  of  inputs  can  be  seen  in  the  DeLaurier  paper  (5); 
however,  in  order  to  compare  the  BEMT  code  to  the  original  DeLaurier  BET  code,  these  same 
inputs  were  used  in  the  BEMT  analysis.  A  comparison  of  the  two  is  shown  in  figure  12. 


BEMT  Thrust  Calculation  vs.  DeLaurier  BET 


Dynamic  Twist 


Figure  12.  BEMT  thrust  vs.  DeLaurier  BET  using  DeLaurier  paper  (5)  inputs. 

It  is  seen  that,  for  the  same  set  of  inputs,  the  DeLaurier  BET  model  overestimates  thrust  for 
dynamic  twist  values  of  0  <P<  77m,  while  the  BEMT  model  calculates  higher  values  of  thrust 
for  a  dynamic  twist  of  [3  >  77m.  It  is  difficult  to  draw  direct  conclusions  from  the  specific  values 
of  thrust  in  this  chart,  since  the  input  [3  is  not  a  real  geometric  input  (e.g.,  not  span,  chord,  etc.) 
and  is  only  used  to  estimate  the  amount  of  wing  twist.  However,  it  can  be  seen  that 
implementation  of  BEMT  still  follows  the  same  trends  when  using  the  same  set  of  inputs  as  the 
DeLaurier  model,  which  verifies  that  the  BEMT  model  is  behaving  “correctly”  (assuming  the 
original  DeLaurier  model  was  correctly  validated). 

The  above  plot  of  BEMT-calculated  thrust  versus  dynamic  twist  shows  the  thrust  values 
converged  upon  between  momentum  theory  and  the  DeLaurier  BET  model  in  the  loop  shown  in 
figure  13.  These  converged  values  come  from  the  minimum  error  between  BET-calculated  thrust 
and  estimated  thrust  used  in  momentum  theory.  Convergence  of  this  BEMT  loop  was  treated  as  a 
minimization  or  optimization  problem,  where  error  is  a  function  that  must  be  minimized. 
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Error  —  ( Calculated  Thrust)BEr  —  Estimated  Thrust 


(5) 


Many  optimization  algorithms  exist;  however,  in  order  to  save  time  on  implementing  a  more 
advanced  optimization  routine,  the  approach  in  this  project  was  simply  to  use  a  direct  search 
method.  In  direct  search  method  algorithms,  the  function  is  evaluated  using  multiple  inputs  over 
a  specified  range,  and  the  resulting  minimum-value  output  is  known  to  be  the  optimum.  A  plot  of 
convergence  upon  the  thrust  value  in  the  BEMT  loop  is  shown  in  figure  13.  It  is  seen  that  the 
error  function  is  nonlinear  and  contains  multiple  local  minima;  when  searching  for  the  global 
minimum  of  a  function,  multiple  local  minima  can  cause  the  optimization  algorithm  to  converge 
on  a  local  minima  should  the  incorrect  input  range  be  used  (function  inputs  bounded  around  a 
local  minima  as  opposed  to  the  global  minimum). 


Figure  13.  Convergence  of  thrust  in  BEMT  loop. 

In  practical  terms,  this  means  that  the  BEMT  analysis  could  converge  upon  the  wrong  thrust 
estimate.  Therefore,  it  is  currently  necessary  for  the  BEMT  code  user  to  exercise  discretion  when 
deciding  upon  the  range  of  inputs  for  thrust  estimation.  This  requires  the  user  to  have  an  idea  of 
where  the  converged  value  should  lie,  e.g.,  experimental  data  can  be  used  in  conjunction  with 
other  methods  of  calculating  thrust  to  search  for  a  minimum  near  the  correct  location.  In  the 
future,  a  more  robust  optimization  routine  should  be  used  that  doesn’t  require  the  user  to  employ 
intuition  about  the  “correct”  range  of  inputs  and  outputs. 

One  main  component  of  the  project  was  to  use  3-D  wing  scatter  plots  to  calculate  the  lift 
coefficient  of  various  wing  airfoil  slices  along  the  span  throughout  the  flapping  cycle.  In  the 
BEMT  code,  airfoil  shapes  are  analyzed  at  discrete  points  along  the  span,  and  thin  airfoil  theory 
is  used  to  estimate  Q  for  the  airfoil.  Airfoil  slices  taken  from  the  3-D  wing  scatter  plots  at  the 
root,  mid-chord,  and  tip  are  shown  in  figure  14.  The  lines  are  a  3  order  polynomial 
approximation  of  the  scatter  coordinates,  using  the  procedure  shown  in  figure  1 1 .  Note  that  the 
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magnitude  of  curvature  shown  is  a  slight  exaggeration  due  to  the  step  size  on  the  z-axis  (max 
thickness  is  much  smaller  than  shown).  It  is  seen  that  the  wing  will  be  deformed  both  concave  up 
and  concave  down  along  the  span,  and  these  deformations  will  change  throughout  the  flapping 
cycle. 


Figure  14.  Wing  airfoil  approximations  across  the  span. 

3.2  3-D  Wing  Scatter  Plot  Airfoil  Representation 

Using  the  same  airfoil  shapes  seen  above  for  calculations  in  thin  airfoil  theory  (figure  15),  it  is 
seen  that  the  Q  versus  a  curve  varies  greatly  for  each  of  the  airfoil  shapes.  Generally,  more 
curvature  when  concave  down  will  result  in  more  positive  values  of  Ci,  as  well  as  a  higher  value 
of  zero  lift  angle  of  attack.  This  is  reflected  in  figure  15,  where  the  tip  shape  is  the  most  concave 
down  and  produces  the  highest  values  of  Q,  while  the  mid-chord  shape  (concave  up)  produces 
negative  Q  values  at  most  angles  of  attack.  It  is  interesting  to  note  the  predicted  Ci  values  at  high 
angles  of  attack  (AoA)  (e.g.,  15+),  where  the  highest  estimated  value  for  the  three  airfoils  shown 
is  Ci  ~3.4  at  AoA  =  20°  (tip  airfoil),  and  even  the  lowest  predicted  Ci  at  20°  is  for  the  mid-chord 
airfoil,  where  Q  ~  0.6.  Because  the  RoboRaven  cruises  at  AoA  =  20°,  it  is  likely  that  these  high 
AoA  Ci  calculations  using  thin  airfoil  theory  are  overestimates  (this  is  expanded  upon  later). 
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Figure  15.  Q  vs.  a  curves  across  span. 

3.3  BEMT  Code  Lift  and  Thrust  Calculations 

As  lift  and  thrust  are  the  two  main  forces  analyzed  in  this  project,  one  interesting  comparison  of 
the  entire  BEMT  code  (including  wing  scatter  data)  is  the  change  in  both  thrust  and  lift  versus 
AoA.  This  is  shown  in  figure  16,  When  calculating  thrust  and  lift  versus  AoA,  the  lift  increases 
roughly  linearly  with  AoA.  This  makes  sense  physically  and  analytically,  because  the  lift 
coefficient  was  calculated  using  thin  airfoil  theory.  The  thrust  change  with  AoA  is  nonlinear  and 
appears  to  be  relatively  unpredictable  through  the  range  of  AoAs  shown.  According  to 
momentum  theory,  when  a  wing  is  placed  at  an  increasing  AoA,  a  correspondingly  increasing 
amount  of  airflow  will  be  deflected  downward,  thereby  increasing  lift.  As  this  downward 
deflection  increases  (AoA  increases),  more  airflow  will  be  used  for  lift  and  less  for  thrust 
creation,  thus  thrust  should  generally  decrease.  Conversely,  as  AoA  decreases  into  negative 
values,  the  wing  flapping  should  move  more  air  in  the  thrust  plane.  This  concurs  with  figure  16, 
where  thrust  appears  to  reach  a  maximum  at  approximately  AoA  =  -2°. 

It  is  interesting  to  note  that  the  BEMT  thrust  appears  to  remain  relatively  unchanged  between 
~2°  <  AoA  <  12°.  This  differs  greatly  from  the  linear  decrease  in  thrust  from  AoA  =  12°  to 
AoA  =  25°.  In  comparison  to  the  BEMT  thrust,  the  DeLaurier-calculated  thrust  appears  to  be 
roughly  parabolic,  with  a  maximum  thrust  around  AoA  =  5°.  For  both  the  DeLaurier  and  BEMT 
calculations  shown,  the  set  of  inputs  were  exactly  the  same,  e.g.,  P  =  07m  and  the  same  3-D 
wing  scatter  plots  were  used  for  wing  shape  throughout  flapping. 
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BEMT  Thrust 
BEMT  Lift 
DeLaurier  Lift 
DeLaurier  Thrust 


Figure  16.  Calculated  thrust  and  lift  vs.  AoA. 

For  the  experimental  wind  tunnel  measurements  of  the  RoboRaven,  it  was  found  that  average  lift 
and  thrust  were  2.38  N  and  1.14  N,  respectively,  when  oriented  at  a  20°  AoA  in  the  wind  tunnel. 
This  AoA  was  chosen  because  the  RoboRaven  was  measured  to  cruise  in  free  flight  at  20°  with 
the  horizontal  free  stream.  The  BEMT  code  calculates  lift  and  thrust  at  6.44  N  and  -2.62  N, 
respectively,  for  an  input  of  AoA  =  20°  (table  1).  It  is  likely  that  the  lift  is  overestimated 
significantly  because  of  the  underlying  limitations  of  thin  airfoil  theory.  Thin  airfoil  theory 
calculates  Q  as  entirely  linear  regardless  of  AoA,  whereas  the  stall  region  at  high  angles  results 
in  a  sudden  drop  in  Q. 


Table  1.  Lift  and  thrust  for  (5  =  0,  a  =  20°. 


Experimental 

Measurements 

BEMT  Calculations 

DeLaurier  BET 

Average  lift  [N] 

2.38 

6.44 

6.90 

Average  thrust  [N] 

1.14 

-2.62 

-1.50 

Thrust  and  Lift  vs.  AoA 

in 

Q 

c 

z  4 

flj  9 

O  ,  ,  X 

LJ-  u 

-10  -5 

~  - 1 - 1 - 

5  10  15  20  25  3 

A 

CL 

Angle  of  Attack  [Degrees] 

Comparing  these  experimental  measurements  with  figure  16  showing  BEMT  calculated  thrust 
and  lift  over  a  range  of  AoAs,  the  closest  match  between  these  forces  is  at  AoA  =  10°,  where 
lift  =  2.51  N  and  thrust  =  0.20  N.  Lift  calculation  via  thin  airfoil  theory  appears  to  be  an 
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overestimate  at  high  AoAs,  since  the  actual  RoboRaven  wind  tunnel  measurements  were 
performed  at  AoA  =  20°  and  match  the  BEMT  thrust  calculations  at  AoA  =  10°.  Additionally, 
the  BEMT  code  significantly  underestimates  thrust,  especially  at  high  AoAs.  The  thrust 
estimation  at  AoA  =  20°  is  much  lower  than  the  experimental  results,  especially  when  comparing 
the  BEMT  results  to  the  original  DeLaurier  BET  results.  It  is  seen  that  the  DeLaurier  BET  is 
more  positive  at  AoA  =  20°,  but  is  still  very  different  from  the  experimental  measurement 
(Thrust  =  —1.5  N  for  DeLaurier  vs.  1.14  N  experimentally).  It  is  clear  that  thrust  effects  are  not 
being  accounted  for  in  both  DeLaurier  and  the  BEMT  model,  which  are  causing  both  to 
underestimate  thrust  at  higher  AoAs. 


4.  Summary  and  Conclusions 


The  initial  goals  of  this  project  were  to  improve  upon  existing  aerodynamic  models  for 
FWUAVs  by  implementing  BEMT  and  more  accurately  depict  wing  shape  in  flight  by  using 
experimental  wing  deformation  measurements.  The  first  goal  was  accomplished  by  combining 
J.D.  DeLaurier’s  BET  model  with  Shkarayev’s  momentum  theory  model  and  iterating  on 
induced  velocity  to  converge  on  thrust  between  the  two  models.  The  second  goal  was 
accomplished  through  the  use  of  experimentally  obtained  3-D  scatter  plots  of  the  FWUAV  wing 
shape  throughout  the  flapping  cycle.  Many  observations  and  conclusions  can  be  drawn  from  the 
general  methodology  used  in  this  project,  as  well  as  on  the  specific  method  used  for 
implementation. 

The  main  reason  for  implementing  a  BEMT  analysis  is  that  convergence  between  two  different 
calculation  methods  and  using  momentum  theory  for  induced  velocity  calculation  will  yield  a 
more  accurate  final  result.  However,  scrutiny  must  still  be  given  to  the  original  BET  and 
momentum  theory  models  that  were  combined  to  create  the  BEMT  model,  as  each  model  will 
still  have  individual  limitations.  In  this  project,  the  DeLaurier  BET  model  approximates  wing 
flapping  motion  as  generally  sinusoidal,  with  inputs  such  as  flapping  angle,  flapping  frequency, 
and  wing  geometric.  While  this  general  method  does  resemble  actual  flapping  in  its  periodic 
motion,  the  3-D  scatter  plots  from  DIC  data  used  in  this  project  can  be  used  to  more  accurately 
represent  wing  flapping  shapes.  Because  wing  shape  is  very  important  when  calculating  lift  and 
thrust,  further  improvements  to  wing  flapping  motion  modeling  should  be  made.  DIC  can  be 
used  to  measure  wing  flapping  velocity  and  acceleration;  creating  function  regressions  from 
these  DIC  measurements  could  be  used  in  place  of  purely  sinusoidal  motion.  Ideally,  this  would 
enable  the  aerodynamic  model  to  have  a  coupling  between  both  wing  deformation  and  motion  to 
be  used  in  BEMT  calculations. 
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Similarly,  the  Shakarayev  model  used  in  this  analysis  was  relatively  simple,  as  described  in 
section  3,  since  flow  properties  are  assumed  to  be  uniform  throughout  the  entire  disk  area. 
Because  flight  forces  across  the  wing  vary  greatly  with  wing  shape  and  speed  (tip  speeds  are 
higher  than  root  speeds),  an  improved  momentum  model  would  account  for  these  differences 
across  the  span.  This  approach  has  been  used  for  rotorcraft,  where  both  models  are  functions  of 
geometric  properties,  as  well  as  spanwise  location  along  the  rotor  blades. 

Another  main  component  of  the  project  was  to  use  thin  airfoil  theory  to  approximate  the  lift 
coefficient  across  the  wing  throughout  the  flap  cycle.  While  distinct  in  its  approach  (the  original 
DeLaurier  BET  model  does  not  use  experimental  wing  shape  measurements),  the  resulting  lift 
coefficient  calculations  seem  to  largely  overestimate  the  true  wind  tunnel  load  cell 
measurements.  This  is  attributed  to  the  ideal  nature  of  thin  airfoil  theory — thin  airfoil  theory 
calculates  a  purely  linear  Ci  versus  a  curve,  which  means  the  user  could  technically  input  an 
extremely  large,  unrealistic  AoA  (e.g.,  AoA  =  50°)  and  thin  airfoil  theory  would  calculate  a 
correspondingly  high  Q  with  no  stall  effects  included.  Thus,  while  the  usage  of  DIC  wing  shape 
measurements  is  one  step  toward  a  better  approximation  of  lift  coefficient  for  a  flapping  wing,  it 
would  be  beneficial  to  improve  the  high  AoA  Q  calculations.  While  it  is  possible  that  the  actual 
lift  coefficient  at  high  AoAs  for  the  most  concave  portions  of  the  wing  truly  are  in  the 
neighborhood  of  Q  =  3-4,  it  is  unlikely.  Better  accounting  for  separation  effects  at  high  AoA 
(where  many  FWUAVs  fly)  should  be  a  focus  in  the  future. 

The  final  goal  of  this  project  was  to  improve  the  thrust  estimation  within  the  BEMT  code.  When 
comparing  the  BEMT  results  with  the  experimental  wind  tunnel  results  at  AoA  =  20°,  it  was 
found  that  the  BEMT  code  significantly  underestimates  thrust  at  higher  AoAs.  The  exact  reason 
why  this  is  occurring  has  not  yet  been  defined;  however,  the  DeLaurier  BET  code’s  function 
approximation  of  thrust  (which  the  BEMT  code  is  based  on)  appears  to  underestimate  thrust  and 
overestimate  lift.  It  is  believed  the  code  should  be  modified  for  these  issues  (this  project  does  not 
explore  how  to  fix  the  problems),  such  that  the  forces  on  an  airfoil  section  better  estimate  the 
thrust  generated.  Since  the  BEMT  code  was  created  as  an  attempt  to  improve  thrust  calculation, 
it  appears  the  specific  BEMT  method  used  in  this  project  should  be  improved  further  since  the 
thrust  is  underestimated  more  than  when  using  the  original  DeLaurier  BET  code.  One  possible 
cause  is  that  the  Shkarayev  momentum  theory  approximates  the  partial  disk  area  created  by  wing 
flapping  as  having  uniform  flow  properties.  In  reality,  since  the  flow  changes  significantly  along 
the  wing  span,  these  changes  need  to  be  accounted  for  in  momentum  theory.  This  type  of 
momentum  theory,  which  is  based  on  the  radial/spanwise  location,  has  been  used  in  rotorcraft 
theory  for  decades  (see  Leishman  [S]),  and  an  adaptation  for  flapping  wings  could  prove 
beneficial. 


16 


5.  References 


1.  Aerovironment,  Inc.  www.avinc.com  (accessed  31  July  2013). 

2.  Young  Australian  at  the  Forefront  of  UAV  Innovation.  (2009,  August).  Retrieved  2013,  from 
Australian  Government  Defense  Materiel  Organisation: 

http://www.defence.gov.au/dmo/news/ontarget/2009/iul09/inn.cfm  (accessed  31  July  2013). 

3.  Lehman,  F.;  Pick,  S.  The  Aerodynamic  Benefit  of  Wing-Wing  Interaction  Depends  on 
Stroke  Trajectory  in  Flapping  Insect  Wings.  J.  of  Experimental  Biology  2007,  210. 
http://ieb.biologists.Org/content/2 10/8/1 362.full.pdf +html. 

4.  University  of  Maryland,  http://youtu.be/miOWpwbnmTw  (accessed  31  July  2013). 

5.  DeLaurier,  J.  D.  An  Aerodynamic  Model  for  Flapping-Wing  Flight.  Aeronautical  Journal 
1993,  Institute  of  Aerospace  Studies,  Ontario. 

6.  Shkarayev,  S.;  Silin,  D.  Applications  of  Actuator  Disk  Theory  to  Membrane  Flapping 
Wings.  AIAA  Journal  October  2010,  48  (10). 

7.  Anderson,  J.  Fundamentals  of  Aerodynamics',  McGraw-Hill  Inc.:  New  York,  1984. 

8.  Leishman,  J.  G.  Principles  of  Helicopter  Aerodynamics.  Cambridge  University  Press.  2006. 


17 


1  DEFENSE  TECHNICAL 

(PDF)  INFORMATION  CTR 
DTIC  OCA 

1  GOVT  PRINTG  OFC 
(PDF)  A  MALHOTRA 

2  DIRECTOR 

(PDF)  US  ARMY  RESEARCH  LAB 
RDRL  CIO  LL 

IMAL  HRA  MAIL  &  RECORDS  MGMT 
1  DIRECTOR 

(PDF)  US  ARMY  RESEARCH  LAB 
RDRL  VTV  J  GERDES 


18 


